#!/usr/bin/python
import sys
import os
import matplotlib.pyplot as plt
import numpy as np

# parameters.
file1 = sys.argv[1]
file2 = sys.argv[2]
wsize = int(sys.argv[3])
gap = int(sys.argv[4])


def read_window(fpath, windsz):
	# read data.
	fin = open(file1, "rb")
	lines = fin.readlines()
	fin.close()
	
	# initalize coverage array.
	y = np.zeros(len(lines), dtype=np.float)
	
	# calulate average coverage.
	avg_list = []
	i = 0
	for line in lines:
		if line[0] == "#": continue
		
		# get coverage.
		tmp = line.split()

		# get single and double.
		single = float(tmp[1])
		single += float(tmp[2])
		single += float(tmp[4])
		single += float(tmp[5])
		
		double = float(tmp[6])
		double += float(tmp[7])
			
		# get ratio.
		if double == 0.0:
			ratio = 0.0
		else:
			ratio = single/double
			
		# add to list.
		if len(avg_list) < windsz:
			# initialize window.
			avg_list.append(ratio)
			y[i] = 0.0
		else:
			# put current, pop last.
			avg_list.pop()
			avg_list.append(ratio)
			
			# calc average.
			tot = 0.0
			for t in avg_list:
				tot += t
				
			#print "%i\t%f" % (i, tot/len(avg_list))
			y[i] = tot/len(avg_list)
		
		i += 1
	return y
	
# get coverages.
y1 = read_window(file1, 25)
y2 = read_window(file2, 25)
ygap = np.zeros(gap, dtype=np.float)
for i in range(len(ygap)):
	ygap[i] = 1.0

# make x-axis.
x1 = np.zeros(len(y1), dtype=np.int)
x2 = np.zeros(len(y2), dtype=np.int)
xgap = np.zeros(len(ygap), dtype=np.int)

# fill x axis appropriate.
j = 0
k = 0
for i in range(len(x1)):
	x1[j] = k
	j += 1
	k += 1
j = 0
for i in range(len(xgap)):
	xgap[j] = k
	j += 1
	k += 1
j = 0
for i in range(len(x2)):
	x2[j] = k
	j += 1
	k += 1
	
	
# plot the graph.
plt.plot(x1, y1, "r-", xgap, ygap, "bs", x2, y2, "g^")
plt.ylabel("single cov/double cov")
plt.show()
